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ABSTRACT 

We investigate the orbital stability of a putative Jovian planet in a compact binary v Oc- 
tantis reported by Ramm et al. We re-analyzed published radial velocity data in terms of 
self-consistent Newtonian model and we found stable best-fit solutions that obey observa- 
tional constraints. They correspond to retrograde orbits, in accord with an earlier hypothesis 
of Eberle & Cuntz, with apsidal lines anti-aligned with the apses of the binary. The best-fit 
solutions are confined to tiny stable regions of the phase space. These regions have a structure 
of the Arnold web formed by overlapping low-order mean motion resonances and their sub- 
resonances. The presence of a real planet is still questionable, because its formation would be 
hindered by strong dynamical perturbations. Our numerical study makes use of a new com- 
putational Message Passing Interface (MPI) framework Mechanic developed to run massive 
numerical experiments on CPU clusters. 
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1 INTRODUCTION 



The v Octantis is a single-line spectroscopic binary composed of 
the v Oct A Kl III giant primary (1.4 + 0.3 M©) and an unseen red 
dwarf secondary v Oct B, K7-M1 V (0.5 + 0.1 M Q ) separated by 
~ 2.55 + 0.13 au. Colacevich (1935 ) collected 11 radial velocities 
(RVs) of the primary, and Alden ( 1939 ) determined the astrometric 
orbit. Recently, the orbital period P b in = 1050.11 + 0.13 days and 
eccentricity e h = 0.2358 + 0.0003 were determined by Ramm ^TaT] 
(2009 ). Remarkably, they also derived the inclination / bin = 71° 
with an error less than 1° on the basis of Hipparchos astrometry and 
new 223 precision RVs. Ram m et al.| ( |2009| > discovered residual 
variability of these RVs at a level of ~ 200 ms _1 and explained 
this by an S-type (Dvorak 1984 ) Jovian planet, v Oct Ab, having 
a minimal mass m p sin i p = 2.5 m Jup in an orbit with a semi-major 
axis a p = 1.2 + 0.1 au and an eccentricity e p = 0.123 + 0.037. 

The putative planet is unusual, because the derived semi-major 
axis implies its orbit roughly in the middle between massive pri- 
mary and secondary. A formation of such a planet in the com- 
pact binary can be hardly explained (e.g., Kley 2010). Analytic 
or general stability criteria formulated for the restricted and gen- 
eral three-body problem like the Hill stability criterion, the results 
of Holman & Wiegert| ft999) , or the resonance-overlap criterion 
by Wisdom (1983) are violated if the planet is assumed to be in a 
prograde orbit. According with the determination of stability limits 
in binaries (Holman & Wiegert 1999), a coplanar, prograde plane- 
tary configuration might be stable only up to astrocentric distance 
~ 0.64 au ( |Ramm et al.1[2009] ). Indeed, |Eberle & Cuntz] ( [2010} 
found that the Keplerian best-fit model in the discovery paper is 
strongly unstable in just 1-10 binary periods time-scale for the 
case of prograde motion. This contradicts the planetary hypothe- 



sis, although the discovery paper basically excludes non-planetary 
sources of the observed signal, like stellar spots or pulsations. 

Many explanations of the observed RVs variability are still 
possible, which might resolve the apparent paradox of the strongly 
unstable system. The residual RVs may be implied by stellar chro- 
mo spheric activity, a different number of planetary companions, 
systematic errors of the observations, instrumental instabilities and 
non-Gaussian uncertainties, such as the red noise (Baluev 2011), 
or a different orbital model than the configuration in the discov- 
ery paper. Indeed, Morais & Correia ( 2012 ) consider a hierarchi- 
cal triple system with an unseen binary in the place of secondary 
y Oct B. Such an orbital setup forces a precession of the primary 
orbit around the barycentre of the inner binary and it could mimic 
the RV variations attributed to the putative planet. 

Earlier, |Eberle & Cuntz| ( |2010| ) following the stability stud- 
ies of the restricted three body problem by Jefferys (1974 ) found 
that a retrograde orbit of the planet lies in much wider stable zone 
than the direct configuration. In a new work, |Quarles et ah] ([201 2) 
investigate this model, extending the set of orbital configurations, 
through a numerical study in the framework of the elliptic, re- 
stricted three-body problem. 

In more than a decade, modeling the precision RVs teaches 
us that the dynamical analysis of the v Octantis system done in 
the discovery paper (Ramm et al. 2009), and following studies 
( [Eberle & Cuntz|2010||Quarles et al.|2012| > do not seem fully con- 
sistent with its dynamical character. Because the secondary mass 
may be almost a half of the primary mass, the relatively wide plan- 
etary orbit is strongly perturbed. The perturbation parameter, ex- 
pressed by the mass ratio of the binary, can be as large as ~0.38 
(Ramm et al. 2009). In such a case the mutual interactions are im- 
portant to compute the RV signal (Laughlin & Chambers 2001). 
Furthermore, it was assumed that the whole system is co-planar, 
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no matter if the planetary orbit is prograde or retrograde. However, 
non-planar orbits in compact binaries might appear due to violent 
post-formation scenarios due to the planet-planet scattering (e.g., 
|Adams & Laug hlin 2003). This may imply highly inclined con- 
figurations and their dynamics much complex than in a coplanar 
case (e.g., Migaszewski & Gozdziewski 201 1). The dynamical sim- 
ulatio ns done by|Eberle & Cuntz| ([2010) concern initially aligned 
orbits. |Quarles et al.| ( |2012| ) study a system with initial planetary 
eccentricity e p = and two initial longitudes A p = 0°,180°, re- 
spectively (corresponding to 3 o'clock and 9 o'clock positions in 
their terminology). While a particular initial orientation of the or- 
bits may be protecting factor for the stability, the initial orbital 
phases should not be fixed or selected a 'priori, if the tested config- 
uration is meant to be consistent with the observations, where we 
aim to analyze the dynamics of a real system. As a simple exam- 
ple, consider Keplerian RV signal computed from the well known 
formulae RV(v) = K[cos(a> + v) + ecos(a>)], where K is the semi- 
amplitude, CD is for the pericenter argument of the orbit, e is the 
eccentricity and v is the true anomaly. Function RV(v) changes its 
sign, if instead of co = 0° one sets co = 180° at the initial time. 
Obviously, a more complex modification of RV(v) is introduced if 
other angles are altered, such as orbital longitudes. In general, the 
initial orbital configuration determines the dynamical evolution and 
the reflex motion of the primary, which should be consistent with 
observations. 

In this work, we aim to verify and improve the kinematic 
(Keplerian) model of the v Oct planetary system by searching for 
the best-fit configurations in terms of self-consistent dynamical, N- 
body model (Laughlin & Chambers 2001 and references therein). 
To resolve the fine structure of the phase space, and to investigate 
the dynamics and long-term stability of the planet, we apply the 
fast-indicator MEGNO (Cinc otta"et al.|20 03) adapted to our new 
multi-CPU computing environment Mechanic. In this paper, we 
demonstrate a non-trivial application of this software. It was an- 
nounced by Sl onina et al.|p012| . The Mechanic code is described 
in an accompanying work to this paper (Slonina et al., in prepara- 
tion), and is freely available on-line (http : //git . ca . umk . pi). 

This paper is structured as follows. After this introduction, 
section 2 is devoted to the dynamical analysis of orbital solution 
proposed in ( [Eberle & Cuntz|2010| ) to possibly global extent. We 
attempt to visualize complex resonant structures in the vicinity 
of presumable retrograde planetary configurations. This part has 
a general character, because we analyze unstable behaviours of S- 
type planets in compact binaries. In the next section 3, we focus on 
an re-analysis of the RV observations and on deriving the best-fit 
parameters of the v Oct system in terms of three distinct models 
of the RV: the Keplerian (kinematic) formulation, Newtonian (dy- 
namic) model, and the Newtonian model with stability constraints 
(GAMP, Gozdziewski et al. 2008). In section 4 we demonstrate that 
stable, retrograde orbits consistent with the RVs data may be found, 
indeed, but these three RV models lead to significantly different or- 
bital architectures. Conclusions are given in section 5. 



2 THE DYNAMICS OF THE v Octantis SYSTEM 

To resolve the paradox of the v Octantis system instability, Eberle 
|& Cuntz| ( |2010| ) searched for stable orbits at a grid of primary 
mass (1.1, 1.4 and 1.7) M Q , three mass ratios Qi = 0.2593,0.2754, 
0.2908), and 30-step initial distance ratio p between the planet and 
the secondary, p e (0.22,0.54) = (0.56, 1.38) au. The equations of 
motion were integrated for 10 3 years (~ 350 binary periods P b in) 



for prograde orbits, and for 1 x 10 4 years (~ 3500 P b in) for retro- 
grade orbits. The initial orbits were aligned with both secondaries 
fixed in their apoastrons with respect to the primary. These integra- 
tions confirmed the theoretical stability limit p ~ 0.25 = 0.64 au 
for prograde orbits in (Holman & Wiegert 1999). For the retrograde 
case, the stability limit was found much larg er, indeed, p ~ 0.479 
consistent with the formal 3cr best-fit error in (Ra mm et al.|20 09'p" 

2.1 Stability of a system with initially aligned apsides 

To extend the study of Eberle & Cuntz (2010), we applied the dy- 
namical maps technique. It helps to illustrate the global structure 
of the phase space, as well as to identify possible resonances as 
the origin of strong instability observed in this system. We use 
MEGNO which is a measure of the maximal Lyapunov exponent 
(Cinc otta et al.|2 003) as the stability indicator. Our serial-CPU soft- 
ware to integrate the equations of motions and the variational equa- 
tions of the planetary problem was encapsulated in the Mechanic 
MPI module called the CSP{^] A strongly interacting planetary sys- 
tem with collisional orbits requires an appropriate integrator. We 
applied the Bulirsh-Stoer-Gragg scheme implemented as reliable 
and well tested ODEX routine ( Hai rer et al.|1993|). It pr ovides very 
good accuracy and performance (see also, [Chambers 1999). The 
total integration time of a single initial condition was not less than 
10 4 periods of the binary, and in some cases more than 3 x 10 5 peri- 
ods. This time is typically one-two orders of magnitude longer than 
the time- span of direct numerical integrations in Eberl e & Cuntzj 
(2010), and provides at least 10-100 longer estimate of the La- 
grange stability, in accord corrwith features of MEGNO ( |Cincotta| 
|et al.|2003) . The time scale of MEGNO integrations is long enough 
to detect most significant mean motion resonances ( Gozdziewski 
|et al.|2008| >. 

In the first numerical experiment, we set the primary mass to 
1.4 M© which is equivalent to \i - 0.28. In accord with Ebe rle &| 
|Cuntz| ( [2010| >, the planet and the secondary are placed in the apoas- 
trons of their orbits, hence the initial mean anomalies M p = Ms = 
180°, of the planet and the secondary, respectively. Arguments of 
periastrons and the nodal longitudes were set to 0°. In the follow- 
ing numerical experiments, we varied the inclination of the planet 
(equivalent to the relative inclination of the orbits) between i p = 0° 
(direct coplanar configuration) and i v = 180° (retrograde copla- 
nar configuration). This setup, which we consider as aligned or- 
bits, reproduces the initial configurations in Eberle & Cuntz ( 2010), 
i.e., the 9 o'clock position in their terminolog)^] We computed 
MEGNO maps in high-resolution with a grid of 1440 x 900 pix- 
els that were integrated for 25,000 periods of the binary. The re- 



We found some ambiguity in the calculation of po by |Eberle & Cuntzj 
(2010), who state that in the elliptical case, po denotes the ratio of the ini- 
tial distance of the planet from the primary (a p ) relative to the semi-major 
axis of the binary components {au n )- Following this, for their retrograde 
configuration stable for 10 Myr, we may compute a p (l + e v ) = po«bin, 
where po = 0.379, «bin ~ 2.55 au, and e p = 0.123, hence a p 0.86 au. 
It might be also a p = po«bin - 0.966 au, if the initial apoastron distance 
of the planet was approximated as a p , or, if a p (l + e p ) = po«bin(l + <?bin) 
then a p 1.06 au, which is most likely approximation used in the study of 
jEberle & CunE|(2010} . 

This code is available upon request. 
3 In the literature, the terms of aligned and antialigned periapses (or orbits, 
apsidal lines) are well defined only for direct orbits. Here, for coplanar, 
direct and retrograde orbits these terms are defined through the difference 
of the longitudes of periapses, equal to 0° and 180°, respectively. 
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Figure 1. MEGNO dynamical maps for the geometric setup of the v Octantis planetary system investigated by Eberle & Cuntz (2010). The initial apsidal 
lines of the orbits are aligned. A reference position of the planet in the discovery paper by Ramm et al. is marked with an asterisk. The low-order mean motion 
resonances between the planet and the secondary are labeled. The left panel is for the global view of the phase space for the nominal mass of the secondary, 
the right panel is for an artificial brown dwarf secondary, helpful to identify the MMRs. The raw resolution is 1440 x 900 initial conditions integrated for 
~ 3 x 10 4 periods of the binary. 



suits are illustrated in Fig. [T] in the semi-major axis a v — relative 
inclination j rel -plane for the two secondary masses. The right map 
is for an artificial brown dwarf secondary, to illustrate, through a 
comparison, the strength of the perturbation in the nominal v Oct 
system (the left panel in Fig. [TJ. For prograde configurations, the 
stability limit a v ~ 0.6 au is slightly smaller than its value found 
in the previous papers. For retrograde orbits, it is almost doubled 
and reaches a p ~ 1.2 au. A putative planet could be located at a 
border of this zone, filled with unstable low-order MMRs. A sta- 
ble retrograde orbit found by Eberle & Cuntz (20K)J lies in a re- 
gion spanned by strong 4:1, 3:1, 5:2 and 2:1 MMRs. A close-up 
of this area covering the 5:2 MMR, slightly beyond the semi-major 
axis of the planet determined in (Ramm et al. 2009), is shown in 
Fig. [2] Overlapping of the MMRs and their sub-resonances creates 
a global chaotic zone. This effect is identified as the main source of 
instability of planets in compact binaries (Mudryk & Wu 2006). 



The stability limit depends strongly on the relative inclination 
and the mass of the secondary. Configurations with intermediate 
relative inclinations, roughly between ~ 40° and ~ 140°, and par- 
ticularly solutions close to polar orbits (/ re i ~ 90°) are very chaotic. 
They are associated with the Kozai resonance (Kozai 1962), see 
also section 2.3. For non-coplanar orbits, the border of stable zone 
may be much closer to the primary than 0.4 au. Even for a hypothet- 
ical, low-mass secondary ~14 m Jup (Fig.JT] right panel), polar or- 
bits are very unstable. In this case, stable prograde solutions may be 
found up to ~ 1.2 au. We also note that individual unstable MMRs 
are much wider for prograde orbits than for retrograde configura- 
tions. A theoretical explanation of this phenomenon is given in a 
recent paper by Morais & Giuppone ( 2012). The stability limits are 
governed by the phase- space topology of retrograde and prograde 
MMRs. At the p/q mean motion ratio, the prograde resonance is of 
order p - q while the retrograde resonance is of order p + q, hence 
the resonance order is higher. Although this result is derived for 
the circular, restricted three body problem, it should be valid also 
for small-eccentricity binaries. This is confirmed by Fig. [T]sho wing 
different widths of prograde and retrograde MMRs for £ bin ^0.123. 



2.2 Developing the instability with the secondary mass 

An identification of the origin of instabilities for large mass ratio 
of the binary is not straightforward. Due to overlapping MMRs and 
their sub-resonances, a complex pattern of stable and unstable mo- 
tions appears (compare panels in Fig. [TJ. To study how the unsta- 
ble zones are evolved with increasing perturbation, we performed 
a numerical experiment in which the secondary mass was gradu- 
ally increased from an artificial value of 14 m Jup (a brown dwarf, 
see the right panel in Fig. [TJ, up to the nominal value of 532 m Jup 
(the left panel of Fig.JTJ. For the small mass secondary, the MMRs 
can be easily identified, and we can follow, how they expand when 
the perturbation of the secondary grows. In this way, we may also 
refer to a low-dimensional dynamical system given through model 
Hamiltonian in Froeschle et al. (2000). The dynamics of the S-type 
planet in the binary should reveal qualitatively similar features be- 
cause this is also governed by a perturbed Hamiltonian of the form 
l~i = "Ho + 6*Ki , where "Ho is the integrable (Keplerian) part, and 
9i\ is the perturbation due to mutual interactions of the planet and 
the binary. The perturbation parameter e can be expressed through 
6 ~ m s /m P = pi, where m s is the mass of the secondary companion 
y Oct B and ra P is the mass of the primary v Oct A. 

The results are illustrated through a sequence of dynamical 
maps in Fig. [3] These MEGNO maps were computed for ~ 10 4 
periods of the binary. At each panel, we labeled the lowest-order 
MMRs and the mass of the secondary (perturber). Each Mechanic 
run occupied 256 CPUs, and computations took, depending on ini- 
tial conditions (an extent of the chaotic zone) up to a few hours 
of CPU time. Already for the smallest mass ratio, corresponding 
to brown dwarf mass perturber, a large part of the phase space is 
strongly chaotic (the top-left panel of Fig. [3] also the right panel 
in Fig. [TJ. In the relevant zone of retrograde motions, a pattern of 
narrow MMRs appears. A retrograde orbit might have much more 
stable space to explore than a direct orbit. In this regime, the prob- 
ability of a stable orbit would be close to 1. When the perturber 
mass is increased, the widths of MMRs quickly expand, and for 
m s ~ 380 m Jup (roughly, the left mass limit of the secondary) one 
may observe that strong MMRs are overlapping which emerges a 
zone of global chaos, roughly beyond a v ~ 1 .2 au. For the nom- 
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Figure 2. A close up of the MEGNO dynamical map for the orbital configuration with apsidal lines initially aligned, in a neighborhood of the formal best-fit 
model in (Ramm et al. 2009). The resolution of the left panel is 1440 x 900 pixels and, for a reference, the right panel has the resolution of 200 x 50 pixels. 
Each initial condition was integrated for ~ 3 x 10 4 and ~ 10 4 periods of the binary, respectively. 



inal mass of the perturber (see the bottom right panel in Fig. [3}, 
there are only narrow areas of stable motions. A close-up of that re- 
gion in Fig. [^demonstrates a view of the phase space of the model 
Hamiltonian for large perturbation parameter e (Froe schle et al.| 
2000 ), see also ( Slonina et al. 2012). The fine structures seen in this 
map represent the Arnold web in the three body problem. They ap- 
pear due to overlapping MMRs and their sub-resonances, strongly 
split by the perturbation. This feature is expected because it was 
found in the Outer Solar system ( Guzzo 2005) and, very recently, in 
the Kepler- 11 extrasolar system of six planets (Migaszewski et al. 
|2012). Likely, it has not been shown in the three body problem with 
so much details. The right panel of Fig. [2] shows the MEGNO map 
in the same area, but with lower resolution of 200 x 50 pixels, com- 
puted for 1 x 10 4 periods of the binary. We see that a fine resolution 
1440 x 900 pixels, as well as longer integration time ~ 3 x 10 4 P bin 
are crucial to discover the fine details of the phase space. (See sec- 
tion 4 for more examples of such structures in the neighborhood of 
the Newtonian best-fit model of the v Octantis system, consistent 
with available observations). 

The dynamical stability of planetary orbits in systems exhibit- 
ing strong mutual interactions can be influenced even by very small 
changes of the initial conditions due to the presence of unstable res- 
onances and their overlapping. If the mutual perturbations are small 
enough, and chaotic motions appear on a regular net, they may be 
practically stable over very long times ( Guzzo 2005). Such a state 
of the system is called the Nekhoroshev regime. Otherwise, if the 
chaotic motions do not constitute a regular web, and most of orbits 
form a global chaotic zone, the stability of the system is influenced 
by a strong chaotic diffusion. This regime is related to the reso- 
nance overlapping, and is called the Chirikov regime (Froeschle 
|etal.|2000||Guzzo|2005| >. 

The Arnold web emerging due to three-body and four-body 
MMRs was investigated by |Guzzo| ( [2005 [ |2006| ) in the Outer Solar 
system. He mapped the phase space in the semi-major axes planes 
of Jupiter, Saturn, Uranus and Neptune. Here, we illustrate this fea- 
ture in the three-body problem in the semi-major axis — inclina- 
tion plane, corresponding to a different choice of the canonical ac- 
tions, and appearing due to two-body MMRs. The Arnold web cre- 
ates an intermediate zone between the ordered and strongly chaotic 
motions, which spans relatively wide range of the inclination, ex- 
tending for 10°. In this intermediate region, the phase-space mo- 



tions which are chaotic, may persist over long periods of time. The 
results of this experiment also show that the stability limits for co- 
planar binaries derived numerically by Holman & Wiegert (1999) 
may be strongly affected by the inclination. In our case of moderate 
eccentricity of the planetary orbit, these limits are roughly valid for 
i rel e (0°,40°) and / re i e (140°, 180°). A determination of stabil- 
ity limits in systems with non-zero relative inclinations needs ad- 
ditional extensive numerical study recalling that our computations 
concern a particular initial configuration of the orbits. We postpone 
this subject to a new work. 

2.3 The secular dynamics of the planet 

The origin of a wide strip of unstable motions around polar or- 
bits can be explained with the help of the secular octupole-level 
theory (e.g., Migaszewski & Gozdziewski 2011 and references 
therein). For a direct comparison with the MEGNO integrations, 
we focus again on the initial orbital setup in (Eberl e~& Cuntzj 
|2010| >. For small a p ~ 0.6 au and less, the expansion parameter 
a = a p /a b in ~ 0.25 where the system may be considered as non- 
resonant. This permit us to study the long-term dynamics through 
the secular theory. The range of a p e [0.4, 1.2] au is investigated to 
compare both relevant intervals of direct and retrograde motion al- 
though a becomes as large as 0.5 at the right border of this interval 
and the assumptions of the secular theory may be questionable. 

The results are illustrated in Fig. [4] in the form of dynami- 
cal maps in the (a p , / re i)-plane. The left column is for the nominal 
mass of y Oct B, and the right column is for an artificial, low-mass 
secondary of 14 m Jup . The top row illustrates the maximal eccen- 
tricity attained by the planet over the secular orbital period. Both 
perturbers are massive enough to force extreme eccentricity ~ 1. 
This must imply collisions with the primaries. Indeed, in the mid- 
dle panels showing the maximum apocenter distance of the planet, 
the filled dots mark initial conditions that led to a collision with 
the star. These two zones closely coincide. If the initial semi-major 
axis is increased, the minimum distance of the planet to the sec- 
ondary becomes almost always close to au (see the bottom row 
in Fig. [4]). Filled dots in these panels mark regions, in which angle 
Am = m p - tus librates. 

This experiment reveals that a continuous unstable zone be- 
tween * re i e (40°, 140°) predicted by the secular model closely co- 
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Figure 3. The results of a numerical experiment in which the secondary mass was gradually increased from 14 mj up (the top left panel), through 75 mj up , 
150 mj U p, 250 mj U p, 300 mj up , 380 mj up , 462 mj up , up to the nominal mass 532 mj up of the v OctantisB (the bottom right panel). Initial geometric configuration 
is fixed the same as in the retrograde orbit model by Eberle & Cuntz (2010). A reference position of the planet is marked with an asterisk. Most significant 
MMRs between the planet and the secondary are labeled. The MEGNO integration time is ~ 1 x 10 4 periods of the binary. 



incides with unstable region around the polar orbits exhibited in the 
full three-body system. It means that this unstable region appear 
due to the Kozai resonance which is able to force large eccentrici- 
ties in the secular time-scale. More details and a concise description 
of the Kozai resonance in the context of planetary dynamics may 
be found in a recent review ( Beauge et al. 2012 p. 1065-1068, and 
references therein). Hence the primary source of instability can be 
identified with a secular effect rather than with short-term chaotic 



dynamics associated with the MMRs. Outside unstable zones of the 
Kozai resonance, the primary source of instability are the MMRs 
which appear as discrete areas of chaos for small value of the per- 
turbation parameters, and a zone of global chaos when the pertur- 
bation becomes sufficiently strong and the short-term MMRs may 
govern the global dynamics. Let us note, that although the global 
picture of the secular system is illustrated here for the aligned con- 
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Figure 4. The orbital evolution of the secular system in terms of the octupole-level expansion w.r.t. a = a p /au n . The left column is for the nominal secondary 
me =532 mj up , the right column is for the secondary mass me =14 mj up . The initial geometric elements are fixed to keep Am = m v - ms = 0° in the Laplace 
reference frame (aligned orbits). The top row is for the maximum eccentricity of the planet, the middle row is for the maximum of the apocenter distance of the 
planet and the bottom row is the minimum distance of the planet from the perturber. Dotted areas mark librations of Am in the top panels, collisional orbits 
with the primary in the middle panels, and librations of Am in the bottom panels. 



figuration, it remains qualitatively the same for the anti-aligned or- 
bits. 



3 RE-ANALYSIS OF THE RV DATA 

It is well known that the RV measurements errors originate from the 
formal errors <r m of the Doppler method, systematic uncertainties 
due to instrumental instabilities, as well as from internal chromo- 
spheric variability crj itter of the parent star. The later error is difficult 
to quantify particularly for young or active stars. In the recent lit- 
erature, it is common to derive astrophysical estimates of the jitter 
crjitter ( Wright 2005). Joint error of the measurements are obtained 
by adding both uncertainties in quadrature for each data point. Re- 
cently, [Baluev] ([2009]) has shown how to estimate crj itter as an addi- 
tional/^^ parameter of the fit model. Here, we follow the earlier, 
somehow more traditional approach. In this work, we computed the 
best-fit models after rescaling the measurement errors in (iRamm 



let al.|2009) by two values of ^jitter- For the first class of models, we 
choose crjitter = 5 ms -1 , which is relatively small value, exhibited 
by quiet, non-active stars, and likely too conservative in our case. 
The second family of solutions was obtained for cr-^ = 20 ms" 1 . 
This value of jitter uncertainty seems better justified recalling large 
residuals of the Keplerian model in the discovery paper of a similar 
magnitude (~ 20 ms -1 ). Using these RV data, we optimized differ- 
ent 2-planet models. The astrometric data gathered by Hipparchos 
were used indirectly, by fixing the nodal line and inclination of the 
binary orbit close to the combined astrometric and radial velocity 
solution in (jRamm et al.|2 009). 



3.1 Keplerian (kinematic) model of the RV 

In the first fitting experiment, we search for the best-fit Keplerian, 
non-interacting orbits. Hence, we basically follow the discovery pa- 
per. The best-fit configurations of v Oct system were found after 
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extensive search performed with the hybrid optimization algorithm 
(see Gozdziewski et al. 2008 and references therein). It consists 
on two stages: a quasi-global Genetic Algorithm which provides 
many reasonably precise models, and the Levenberg-Marquardt 
(LM) method which is used to refine these solutions. 

The results are illustrated in Fig. [5] in the form of (xl) 1/2 pa- 
rameter scans illustrated in the orbital period P p - eccentricity e p 
plane of the planet. Curiously, for both selected jitter uncertainties, 
we found a well defined secondary minimum for P p ~ 900 days and 
e p ~ 0.3, besides the best- fit model reported by Ra mm et aL] |2009) 
with P p ~ 420 days and e p ~ 0.1. A choice of the jitter uncertainty 
may lead to significant changes of the (xl) 1/2 shape and alters for- 
mal errors of the best-fit parameters. This effect is important for 
fixing parameter bounds in dynamical studies made after the best- 
fit models are found. Moreover, because the secondary best-fit con- 
figuration has the planetary period comparable with the binary pe- 
riod, its interpretation in terms of the osculating elements becomes 
even more questionable than in the case of the primary minimum 
of (xl) 1/2 > as underlined in the introduction. The secondary solu- 
tion might correspond to 1 : 1 MMR which obviously would require 
a more general Af-body model of the RV which takes into account 
all mutual interactions between the bodies in the system (Laughlin 
|& Chambers|200T] ). 



3.2 The self-consistent, Newtonian model of the RV 

As we have shown above, the Keplerian model of v Oct data is non- 
unique, permitting two different solutions. It is also well known that 
it is independent on orbital inclinations and true masses. Due to the 
massive secondary, the RV-model including mutual interactions be- 
tween all system components might constrain these parameters, in 
spite of a narrow observational window, like in other interacting 
extrasolar planetary systems (Laughlin & Chambers 2001). Hence, 
in the second fitting experiment, to compute the RVs of the pri- 
mary, we integrated the Af-body equations of motion. In this model, 
the primary mass m b i n , its semi-major axis a b in, eccentricity e hin of 
the binary, pericenter argument <^ b in, inclination j bin and the initial 
mean anomaly Mbm of the binary were allowed to vary within the 
formal 3<x error bounds, in accord with the combined RV and as- 
trometric solution in Ra mm et al.|p009| >. The nodal longitude of 
the binary was fixed, £2 bin = 87°. The formal error of this element 
is irrelevant for the dynamics of the system, because selecting the 
nodal line of the binary corresponds to a choice of the x-axis of the 
reference frame. The planet mass m p , its inclination z' p , the nodal 
longitude Q p and the remaining orbital elements with the RV off- 
set were added to the set of 14 free parameters of the fit model. 
The uncertainties of the measurements in ( Ramm et al. 2009) were 
rescaled in quadrature, as in the Keplerian model, to account for the 
intrinsic stellar RV variability. 

Similar to the Keplerian model, we explored the parameter 
space with the hybrid algorithm. After an extensive search, we se- 
lected the best-fit Af-body configuration. Then we reconstructed 
the shape of (xl) 1/2 around this nominal solution in 2-dimensional 
planes of various orbital elements. Its parameters derived with 
^ jitter - 5 ms _1 are given in Table 1 (Fit I). The RV data in 



Ramm 



|eTaT] ( [2009l > consist of many "clumps" of measurements done dur- 
ing one night, which cover less than 0.1 day. Such data points were 
binned, and the resulting data set consists of 91 measurements. This 
step was helpful to improve the CPU-efficiency, due to necessity 
of computing numerical derivatives in the LM algorithm. We ex- 
amined that (^) 1/2 -scans can be reproduced with the full, original 
data set. Obtaining a smooth shapes of (xl) l/2 for non-binned data 



is difficult, due to numerical errors introduced by the numerical 
derivatives. 

The (fl p , £ p )-plane (Fig. [6] the top-lef panel) makes it possible 
to compare the topology of (xl) 112 in both RV-models. Clearly, the 
best-fit Newtonian configuration is qualitatively different from the 
Keplerian models found in the previous section which confirms our 
earlier predictions. The planet eccentricity may be as large as ~ 
0.45. The low estimate of crj itter ~ 5 ms -1 implies two, statistically 
distinct close minima of (xl) l/2 > for a p ~ 1.25 au, and better model 
for a p ~ 1.4 au. These minima are resolved at the formal 2cr-error 
level of (x 2 ) 112 - If tne jitter estimate is larger, likely much more 
realistic, only one solution persist (see Fig. [7] the top left panel). In 
this case, the formal lcr error of a p spans almost 0.3 au. 

Even more interesting results are illustrated in the next panels 
of Figs. [6] and [7] showing scans of ix 2 ) 112 in the semi-major axis- 
inclination, — longitude of the node and — the pericenter argument 
planes, respectively. (xl) l/2 scans computed for <j } 



jitter 



5 ms re- 



veal two minima seen in Fig. [6] Such scans for (7j itter ~ 20 ms -1 
(Fig. [7} illustrate that the global minimum around (a v = 1.40 au, 
e p = 0.45) is well defined in all planes of angular elements. This 
means that even a narrow observational window, covering roughly 
2 planetary periods, makes it possible to estimate all angles of 
the orbit due to strong mutual interactions in the system. This ex- 
periment is a convincing illustration that the Keplerian, kinematic 
model of v Oct is questionable because it looses information of 
the node and inclination of the orbit. A proper estimate of the jit- 
ter uncertainty is vital to derive true errors of the best-fit model. 
We found that the best-fit model (Fit I in Table 1) tends towards 
the lower bound of inclination (69°), with simultaneous increase 
of the secondary mass. This effect appears because the inclination 
of the binary cannot be constrained by the RV data alone. Likely, 
dynamically consistent Af-body model of the v Oct system should 
also incorporate the Hipparchos astrometry, which would be help- 
ful to fix the orbit of the binary fully consistent with all data and 
the three-body model. We postpone this analysis to another work. 

3.3 Newtonian model with stability constraints 

Unfortunately, the best-fit Newtonian model (Fit I in Table 1) found 
through the hybrid search leads to strongly unstable system. Due 
to relatively large parameter error, stable configurations might be 
still possible in its neighborhood. Hence, we conducted an addi- 
tional, extensive search for stable models imposing stability con- 
straints (GAMP algorithm, Gozdziewski et al. 2008). In this ap- 
proach, unstable models are "penalized" by an artificial, large value 
of (x 2 v ) 112 (or an rms) and are excluded during the optimization. The 
results are shown in Fig. [6] Stable fits exhibit systematically larger 
(xl) l/2 than the best-fit model in Table 1. To illustrate the statis- 
tics of such configurations, we over-plotted solutions with an rms 
< 28 ms -1 (white, filled circles) at the top of selected parameter 
scans of (xl) l/2 (Fig- [6}- Let us note that stable models have the 
MEGNO signature ~ 2 for 2000 P hin . This relatively short integra- 
tion period requires an acceptable CPU-overhead from the hybrid 
algorithm, still providing at least 10-100 longer Lagrange stability 
time. Such solutions should survive for 10 4 -10 5 P b in or even longer 
integration time (Cincotta et al. 2003). The statistics in Figure- 
veals that the formal best-fit configurations in terms of "pure" N- 
body model are systematically shifted by more than 3<x from the 
best-fit GAMP (stable) models. This regards all orbital parameters, 
and in particular, the eccentricity and the argument of pericenter. 
This large discrepancy may indicate that the 2-planet model is in 
fact not proper, and should be rejected at all or modified. Curi- 
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ously, we found stable solutions only as retrograde orbits, actually 
in accord with the original hypothesis in (Eberle & Cuntz 2010 ) 
and later investigated by Quar les et al.|p012| >. This is illustrated in 
Fig. [8] showing that differences of the longitudes of periastrons and 
of the nodal longitudes are well bounded around 180°, resulting in 
apsidal lines of the orbits which are significantly anti-aligned (the 
left panel). Large relative inclinations of these stable solutions (the 
right panel) mean retrograde orbits in the range of a p ~ around 
1.2 au. 

To compare the quality of these two families of models in 
terms of an rms, the RV residuals of two selected best-fit mod- 
els are displayed in Fig. [9] the first one derived without stability 
constraints (Fit I in Tab. 1, the left panel), and an low rms stable 
model (the right panel, Fit II in Table 1). In both cases, the residu- 
als exhibit a large magnitude of ~ 120-150 ms _1 , comparable with 
the semi-amplitude of the RV signal ~ 100 ms" 1 due to the putative 
planet itself. In this sense, stable solutions might be still acceptable, 
although their rms is significantly larger from the rms of the formal 
best-fit model configuration. A large variability of the RV residuals 
cannot be well explained by the 1 -planet model. This may indicate 
strong stellar activity mimicking planetary signal. 



3.4 Stability of the best-fit GAMP configurations 

To illustrate the dynamical neighborhood of the best-fit GAMP 
models, we selected two examples of low rms configurations. Their 
osculating elements at the epoch of the first observation in ( |Ramm| 
|et al.|2009] > are given in Table 1 (Fits II and III). Many digits are 
quoted, to reproduce elements of these fit possibly exactly, to make 
it possible to reproduce their stable evolution due to complex and 
chaotic neighborhoods. The formal errors may be estimated graphi- 
cally in Fig.|6]or interpreted as uncertainties of the formal Fit I. The 
primary mass is fixed as equal to 1.4 M©. An rms of this solution is 
~ 25 ms" 1 . 



The top-left panel of Fig. [To| illustrates the MEGNO dynam- 
ical map in the (a p , / p )-plane computed for Fit II (see Tab.l). This 
map shows a global view of the phase space. Note that inclination 
z'p is absolute w.r.t. the astrometric frame , in which the binary in- 
clination is bounded around i hin ~ 71°. This figure reveals that the 
phase space is mostly strongly chaotic. The best-fit configuration 
is found in a narrow stability island close to the 7:2 MMR. For a 
reference, other significant MMRs are labeled. A close-up of this 
map shown in the top-right panel of Fig. [TO] reveals an extremely 
complex structure of this regular island. It spans a tiny region cov- 
ering Aa p ~ 0.02 au and A/ p ~ 20°. In this small island, the phase 
space has again a sophisticated, fractal-like structure of the Arnold 
web. 

The integration time used to derive the global dynamical map 

is 1 X 10 4 P b in- 



in the left panel of Fig. 



10 



Subsequent close-ups 
shown in Fig. 10 were integrated for 2 x 10 4 P b in up to 1 x 10 5 P bin 
per each point, depending on the magnification factor. To detect 
weak MMRs and their sub-resonances, such long integration times 
are indispensable. The resolution of these maps is 1440 x 900 pix- 
els. Calculations were performed with the help of the Mechanic 
installed on the UV SGI supercomputer chimera of the Poznari 
Supercomputing Centre (PolandQ 

We computed a similar set of dynamical maps for another sta- 
ble model, with larger eccentricity (Fit III in Table 1). The results 
are illustrated in Fig. [TT] The top left panel illustrates the dynami- 
cal map in the (a p , £ p )-plane, showing a narrow stable island of this 
solution, the rest of plots are for the (a p , z p )-plane and its close-ups. 



4 Supercomputer chimera has 2048 CPU cores, and runs spanning up to 
24 hrs occupied all available hardware resources. We did not observed any 
bottlenecks in the master-worker inter-communication. This nice perfor- 
mance of the Mechanic code was expected thanks to a small overhead of 
the MPI data broadcasts. These computations validated the software and its 
ability of stable long-runs on a few thousands of CPU cores. 
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Figure 6. The shape of (xl) 1 ^ 2 in different orbital planes of the planet computed in terms of the self-consistent A/-body model of the RV. Jitter uncertainty 
is emitter ~ 5 ms _1 . The RV data in (Ramm et al. 2012) are binned (see the text for details). Quality of solutions is color coded. Smooth, coloured curves are 
for formal l,2,3<x confidence levels of the best-fit model marked with an asterisk. Its elements recomputed on the original data set are given in Table 1 (Fit I). 
Smaller, filled circles mark stable solutions found with through the GAMP search, and projected onto planes of orbital osculating elements at the epoch of the 
first observation in (Ramm et al. 2009). An rms of these solutions is less than 28 ms _1 . 
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Figure 7. The shape of (xl) 1/2 in terms of the Newtonian model illustrated in different orbital parameter planes, as derived for jitter uncertainty cr^ 
20 ms _1 . Quality of solutions is color coded. The best-fit configuration is marked with an asterisk. 



To reveal the Arnold web in more detail, we computed a number of 
close-ups shown in panels of Figs.[l0j{T2] These maps reveal many 
weak sub-resonances, which form a regular net present at any mag- 
nification level. These structures are particularly well seen in the 
bottom-right panel of Fig. [12] They closely resemble the Arnold 



web of the model Hamiltonian (Froeschle et al. 2000 ). Likely, it 
has been never illustrated with such a sharpness and details. We 
also tested other solutions from the GAMP-derived set. The results 
are similar. The best-fit models which obey reasonably well the ob- 
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150 160 170 180 190 200 210 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 

£2 p -n bin [deg] semi-major axis a p [au] 

Figure 8. The statistics of stable best-fit models (see Fig. 6) projected onto selected planes of the orbital osculating elements at the epoch of the first observation 
in (Ramm et al. 2009). These fits are gathered in the hybrid GAMP search (see the text for details). An rms of these solutions is less than 28 ms _1 . 



Table 1. The best-fit astro-centric, osculating Keplerian elements of the v Octantis planetary system at the first epoch in (Ramm et al. 2009). The RV data 
errors in the discovery paper are rescaled with the jitter uncertainty <Xji tter = 5 ms _1 in quadrature. The nodal longitude of the binary orbit is fixed to 87 deg. 
Mass of the primary v Oct A is 1.4 m . Parameters of stable solutions (Fits II and III) are quoted with many digits making it possible to demonstrate stability 
of these solutions in very narrow stable zones. Formal errors of these fits can be estimated as equal or similar to the uncertainties quoted for parameters of 
formal, A/-body Fit I. 





Fit I (N-body, unstable) 


Fit II (GAMP, stable) 


Fit III (GAMP, stable) 


Parameter 


planet v OctAb 


v OctB 


planet v OctAb 


y OctB 


planet v OctAb 


vOctB 


mass m [mj up ] 


3.54 ± 0.07 


572 ± 15 


1.92596 


560.62366 


2.135457 


559.2707 


semi-major axis a [AU] 


1.41 ±0.02 


2.533 ± 0.005 


1.16365 


2.52813 


1.160281 


2.526478 


eccentricity e 


0.40 ± 0.01 


0.2331 ±0.0006 


0.13139 


0.23881 


0.247587 


0.238713 


inclination i [deg] 


84.9 ± 2.2 


69.0 ± 3.3 


110.21102 


71.28090 


102.4152 


71.681 


longitude of the node Q [deg] 


289.2 ± 4.2 


87.0 (fixed) 


262.29811 


87.0 (fixed) 


260.0224 


87.0 (fixed) 


argument of pericenter oj [deg] 


340.6 ± 2.6 


75.0 ±0.1 


127.26299 


74.59137 


88.86501 


74.85992 


mean anomaly M(to) [deg] 


289.2 ± 4.2 


339.5 ± 0.1 


133.16327 


339.762973 


165.3685 


339.5231 



0^)1/2 2.79 3.47 3.44 

RV offset V [m s" 1 ] -6417 ± 1 -6473.682 -6475.317 

rms[ms _1 ] 19.9 25.2 24.8 



] 
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1.1530 1.1535 1.1540 1.1545 1.1550 

planetary semi-major axis a p [au] 

Figure 12. A close-up of a stability island of the best-fit, stable GAMP 
Fit II (Table 1) illustrated in Fig.[l0| The resolution is 1440 x 900 pixels. 
The integration time per pixel is 3 x 10 5 orbital periods of the binary. 



servational constrains are confined to only narrow stability islands. 



4 CONCLUSIONS 

In this work, we aim to improve the original Keplerian models of 
the v Oct system. We re-analysed the RV observations in ( |Ramm| 



|et al.|2009] >. The results confirm our predictions on a significance 
of the mutual interactions between the secondary and the putative 
planet. These interactions are so strong that - if the planet is real - it 
might be possible to derive the orientation of the orbit with respect 
to the binary, and its true mass, in spite of a narrow observational 
window spanning less than 2 orbital periods. 

We confirm that stable retrograde solutions in the v Oct plan- 
etary system are possible as proposed in ( [Eberle & Cu ntz 2010, 
Quarles et al. 2012 ). However, the observational constrains imply 
such orbits with apsidal lines initially anti-aligned with the apses 
of the binary. The long-term stable solutions may be found only in 
relatively small islands of regular, quasi-periodic motions exhibit- 
ing a complex structure of the Arnold web. Stable models must be 
significantly shifted in the parameter space from the Af-body and 
Keplerian best-fit models. It means that the 2-planet, S-type con- 
figuration is most likely inadequate to fit the available data set. It 
remains highly uncertain how a massive, Jovian planet could be 
trapped in such tiny stable regions, or how it might form in glob- 
ally unstable dynamical environment. The observational RV signal 
is very noisy, and the best-fit models reveal a large scatter of resid- 
uals having amplitude comparable with the RV signal itself. These 
arguments make the hypothesis of Jovian planet in orbit around 
y Oct A questionable. Actually, this in accord with the doubts raised 
in the discovery paper (Ramm et al. 2009}. New observations are 
required to confirm or withdraw that explanation of the observed 
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Figure 9. Residuals to representative self-consistent, Newtonian best-fit models illustrated in Fig. [8] The left panel is for the formal, mathematical A/-body 
solution (Fit I in Table 1). The right panel is for a stable, low-rms GAMP-derived model (Fit II in Table 1). 
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Figure 10. Dynamical MEGNO maps of a stable model found in the GAMP search illustrated in the (a p , i p )-plane. The top left panel is for the global view 
of the phase space, the top right map is for stability island of the solution. Subsequent maps are for close-ups of the neighborhood of the best-fit solution. 
Osculating, astrocentric elements of the planet and the secondary are given in Table 1 (Fit II). The nominal elements are marked with the star symbol. The raw 
resolution of the maps is 1440x900 data points. 



RV residual signal, and to search for alternate dynamical models of 
the detected variability. 

An efficient analysis of the observations of extrasolar plane- 
tary systems involves not only different observational techniques 
and data sources, but also orbital optimization algorithms, com- 
bined with analytical and numerical stability studies. In this work, 
we propose and apply our new MPI-based task management sys- 
tem Mechanic to perform high-resolution dynamical mapping of 
the phase space. This framework is devoted to numerical analysis 
of large sets of initial conditions, like the long-term integrations of 
the equations of motion. We work on other applications, like quasi- 
global optimization with the Genetic Algorithms of different obser- 



vational models of the extrasolar planetary systems. The Mechanic 
code with detailed technical documentation and sample modules 
is freely available at the project website ( [Slonina & Gozdziewski 

|20T2] ). 
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Figure 11. Dynamical MEGNO maps of a stable model with moderate eccentricity of the planet as found in the GAMP search illustrated in the (a p ,e p ) and 
(a p , / p )-plane. Bottom panels show close-ups of the top right panel. Osculating, astrocentric elements of the planet and the secondary are given in Table 1 (Fit 
III). The nominal elements are marked with the star symbol. The raw resolution of the maps is 1440x900 data points. 
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